home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / balgen / scaleg.f < prev   
Text File  |  1996-07-19  |  7KB  |  237 lines
  1.       subroutine scaleg (n,ma,a,mb,b,low,igh,cscale,cperm,wk)
  2. c
  3. c     *****parameters:
  4.       integer igh,low,ma,mb,n
  5.       double precision a(ma,n),b(mb,n),cperm(n),cscale(n),wk(n,6)
  6. c
  7. c     *****local variables:
  8.       integer i,ir,it,j,jc,kount,nr,nrp2
  9.       double precision alpha,basl,beta,cmax,coef,coef2,coef5,cor,
  10.      *                 ew,ewc,fi,fj,gamma,pgamma,sum,t,ta,tb,tc
  11. c
  12. c     *****fortran functions:
  13.       double precision dabs, dlog10, dsign
  14. c     float
  15. c
  16. c     *****subroutines called:
  17. c     none
  18. c
  19. c     ---------------------------------------------------------------
  20. c
  21. c     *****purpose:
  22. c     scales the matrices a and b in the generalized eigenvalue
  23. c     problem a*x = (lambda)*b*x such that the magnitudes of the
  24. c     elements of the submatrices of a and b (as specified by low
  25. c     and igh) are close to unity in the least squares sense.
  26. c     ref.:  ward, r. c., balancing the generalized eigenvalue
  27. c     problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981,
  28. c     141-152.
  29. c
  30. c     *****parameter description:
  31. c
  32. c     on input:
  33. c
  34. c       ma,mb   integer
  35. c               row dimensions of the arrays containing matrices
  36. c               a and b respectively, as declared in the main calling
  37. c               program dimension statement;
  38. c
  39. c       n       integer
  40. c               order of the matrices a and b;
  41. c
  42. c       a       real(ma,n)
  43. c               contains the a matrix of the generalized eigenproblem
  44. c               defined above;
  45. c
  46. c       b       real(mb,n)
  47. c               contains the b matrix of the generalized eigenproblem
  48. c               defined above;
  49. c
  50. c       low     integer
  51. c               specifies the beginning -1 for the rows and
  52. c               columns of a and b to be scaled;
  53. c
  54. c       igh     integer
  55. c               specifies the ending -1 for the rows and columns
  56. c               of a and b to be scaled;
  57. c
  58. c       cperm   real(n)
  59. c               work array.  only locations low through igh are
  60. c               referenced and altered by this subroutine;
  61. c
  62. c       wk      real(n,6)
  63. c               work array that must contain at least 6*n locations.
  64. c               only locations low through igh, n+low through n+igh,
  65. c               ..., 5*n+low through 5*n+igh are referenced and
  66. c               altered by this subroutine.
  67. c
  68. c     on output:
  69. c
  70. c       a,b     contain the scaled a and b matrices;
  71. c
  72. c       cscale  real(n)
  73. c               contains in its low through igh locations the integer
  74. c               exponents of 2 used for the column scaling factors.
  75. c               the other locations are not referenced;
  76. c
  77. c       wk      contains in its low through igh locations the integer
  78. c               exponents of 2 used for the row scaling factors.
  79. c
  80. c     *****algorithm notes:
  81. c     none.
  82. c
  83. c     *****history:
  84. c     written by r. c. ward.......
  85. c     modified 8/86 by bobby bodenheimer so that if
  86. c       sum = 0 (corresponding to the case where the matrix
  87. c       doesn't need to be scaled) the routine returns.
  88. c
  89. c     ---------------------------------------------------------------
  90. c
  91.       if (low .eq. igh) go to 410
  92.       do 210 i = low,igh
  93.          wk(i,1) = 0.0d0
  94.          wk(i,2) = 0.0d0
  95.          wk(i,3) = 0.0d0
  96.          wk(i,4) = 0.0d0
  97.          wk(i,5) = 0.0d0
  98.          wk(i,6) = 0.0d0
  99.          cscale(i) = 0.0d0
  100.          cperm(i) = 0.0d0
  101.   210 continue
  102. c
  103. c     compute right side vector in resulting linear equations
  104. c
  105.       basl = dlog10(2.0d0)
  106.       do 240 i = low,igh
  107.          do 240 j = low,igh
  108.             tb = b(i,j)
  109.             ta = a(i,j)
  110.             if (ta .eq. 0.0d0) go to 220
  111.             ta = dlog10(dabs(ta)) / basl
  112.   220       continue
  113.             if (tb .eq. 0.0d0) go to 230
  114.             tb = dlog10(dabs(tb)) / basl
  115.   230       continue
  116.             wk(i,5) = wk(i,5) - ta - tb
  117.             wk(j,6) = wk(j,6) - ta - tb
  118.   240 continue
  119.       nr = igh-low+1
  120.       coef = 1.0d0/float(2*nr)
  121.       coef2 = coef*coef
  122.       coef5 = 0.5d0*coef2
  123.       nrp2 = nr+2
  124.       beta = 0.0d0
  125.       it = 1
  126. c
  127. c     start generalized conjugate gradient iteration
  128. c
  129.   250 continue
  130.       ew = 0.0d0
  131.       ewc = 0.0d0
  132.       gamma = 0.0d0
  133.       do 260 i = low,igh
  134.          gamma = gamma + wk(i,5)*wk(i,5) + wk(i,6)*wk(i,6)
  135.          ew = ew + wk(i,5)
  136.          ewc = ewc + wk(i,6)
  137.   260 continue
  138.       gamma = coef*gamma - coef2*(ew**2 + ewc**2)
  139.      +        - coef5*(ew - ewc)**2
  140.       if (it .ne. 1) beta = gamma / pgamma
  141.       t = coef5*(ewc - 3.0d0*ew)
  142.       tc = coef5*(ew - 3.0d0*ewc)
  143.       do 270 i = low,igh
  144.          wk(i,2) = beta*wk(i,2) + coef*wk(i,5) + t
  145.          cperm(i) = beta*cperm(i) + coef*wk(i,6) + tc
  146.   270 continue
  147. c
  148. c     apply matrix to vector
  149. c
  150.       do 300 i = low,igh
  151.          kount = 0
  152.          sum = 0.0d0
  153.          do 290 j = low,igh
  154.             if (a(i,j) .eq. 0.0d0) go to 280
  155.             kount = kount+1
  156.             sum = sum + cperm(j)
  157.   280       continue
  158.             if (b(i,j) .eq. 0.0d0) go to 290
  159.             kount = kount+1
  160.             sum = sum + cperm(j)
  161.   290    continue
  162.          wk(i,3) = float(kount)*wk(i,2) + sum
  163.   300 continue
  164.       do 330 j = low,igh
  165.          kount = 0
  166.          sum = 0.0d0
  167.          do 320 i = low,igh
  168.             if (a(i,j) .eq. 0.0d0) go to 310
  169.             kount = kount+1
  170.             sum = sum + wk(i,2)
  171.   310       continue
  172.             if (b(i,j) .eq. 0.0d0) go to 320
  173.             kount = kount+1
  174.             sum = sum + wk(i,2)
  175.   320    continue
  176.          wk(j,4) = float(kount)*cperm(j) + sum
  177.   330 continue
  178.       sum = 0.0d0
  179.       do 340 i = low,igh
  180.          sum = sum + wk(i,2)*wk(i,3) + cperm(i)*wk(i,4)
  181.   340 continue
  182.       if(sum.eq.0.0d0) return
  183.       alpha = gamma / sum
  184. c
  185. c     determine correction to current iterate
  186. c
  187.       cmax = 0.0d0
  188.       do 350 i = low,igh
  189.          cor = alpha * wk(i,2)
  190.          if (dabs(cor) .gt. cmax) cmax = dabs(cor)
  191.          wk(i,1) = wk(i,1) + cor
  192.          cor = alpha * cperm(i)
  193.          if (dabs(cor) .gt. cmax) cmax = dabs(cor)
  194.          cscale(i) = cscale(i) + cor
  195.   350 continue
  196.       if (cmax .lt. 0.5d0) go to 370
  197.       do 360 i = low,igh
  198.          wk(i,5) = wk(i,5) - alpha*wk(i,3)
  199.          wk(i,6) = wk(i,6) - alpha*wk(i,4)
  200.   360 continue
  201.       pgamma = gamma
  202.       it = it+1
  203.       if (it .le. nrp2) go to 250
  204. c
  205. c     end generalized conjugate gradient iteration
  206. c
  207.   370 continue
  208.       do 380 i = low,igh
  209.          ir = wk(i,1) + dsign(0.5d0,wk(i,1))
  210.          wk(i,1) = ir
  211.          jc = cscale(i) + dsign(0.5d0,cscale(i))
  212.          cscale(i) = jc
  213.   380 continue
  214. c
  215. c     scale a and b
  216. c
  217.       do 400 i = 1,igh
  218.          ir = wk(i,1)
  219.          fi = 2.0d0**ir
  220.          if (i .lt. low) fi = 1.0d0
  221.          do 400 j =low,n
  222.             jc = cscale(j)
  223.             fj = 2.0d0**jc
  224.             if (j .le. igh) go to 390
  225.             if (i .lt. low) go to 400
  226.             fj = 1.0d0
  227.   390       continue
  228.             a(i,j) = a(i,j)*fi*fj
  229.             b(i,j) = b(i,j)*fi*fj
  230.   400 continue
  231.   410 continue
  232.       return
  233. c
  234. c     last line of scaleg
  235. c
  236.       end
  237.